{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-a09a1d0c-1197-4369-851d-f7409bfe896e",
    "output_cleared": false
   },
   "source": [
    "# Phase-space file (PHSP)\n",
    "\n",
    "Use the output of the following simulation:\n",
    "- Folder: linac/\n",
    "- Already computed results in folder output_ref/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "cell_id": "00001-da98f9a5-10a2-4e23-b64b-b7305982a0d4",
    "execution_millis": 297,
    "execution_start": 1605009800829,
    "output_cleared": false,
    "source_hash": "2036bae7"
   },
   "outputs": [],
   "source": [
    "# Tell Jupyter to plot figure right in the page\n",
    "%matplotlib widget\n",
    "\n",
    "# Module with plot capabilities\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Module with scientific computing functions (matrix/vector)\n",
    "import numpy as np                \n",
    "\n",
    "# Module to read root files\n",
    "import uproot\n",
    "\n",
    "# Modules with reading/write folder/file functions\n",
    "import os\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "cell_id": "00002-6af52d4a-45a2-415f-94c9-c9282735e2d9",
    "execution_millis": 3,
    "execution_start": 1605009801128,
    "output_cleared": false,
    "source_hash": "909ee571"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Current Working Directory (CWD) is: \n",
      " /Users/dsarrut/src/gate/dqprm/2021/gate-exercices/linac\n"
     ]
    }
   ],
   "source": [
    "# The following command display the current working directory (where jupyter has been launched)\n",
    "cwd = os.getcwd()\n",
    "print('The Current Working Directory (CWD) is: \\n', cwd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "cell_id": "00003-5822b599-07ee-455c-8cf0-863f2de7f49b",
    "execution_millis": 61,
    "execution_start": 1605009813330,
    "output_cleared": false,
    "scrolled": true,
    "source_hash": "fe0e75b0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Read PHSP object <TTree b'PhaseSpace' at 0x7fc6df7d26d0>\n",
      "PhaseSpace keys:  [b'Ekine', b'Weight', b'X', b'Y', b'Z', b'dX', b'dY', b'dZ', b'TrackID', b'EventID', b'RunID']\n"
     ]
    }
   ],
   "source": [
    "# read a PHSP\n",
    "root_filename = Path('./output_ref/output-PhS-g.root')\n",
    "try:\n",
    "    f = uproot.open(root_filename)\n",
    "except Exception:\n",
    "    print(\"File '\"+root_filename+\"' cannot be opened, not root file ?\")\n",
    "    exit()\n",
    "\n",
    "# Look for a single key named \"PhaseSpace\"\n",
    "k = f.keys()\n",
    "try:\n",
    "    psf = f['PhaseSpace']\n",
    "except Exception:\n",
    "    print(\"This root file is not a PhaseSpace, keys are: \", f.keys())\n",
    "    exit()\n",
    "    psf = f['PhaseSpace']\n",
    "# now, the variable psf contains a root Tree with various branches: Energy, X, Y, Z etc\n",
    "print('Read PHSP object', psf)\n",
    "print(\"PhaseSpace keys: \", psf.keys())\n",
    "# all branches are set in the variable 'a'\n",
    "a = psf.arrays()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "cell_id": "00004-ff22fa1e-88a1-46f4-8fbd-47c761a976e6",
    "execution_millis": 481,
    "execution_start": 1605009815199,
    "output_cleared": false,
    "source_hash": "7756d4f1"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1b7539adecd746af9051a3d1d2ea3552",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc6e0f95760>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot the E\n",
    "nbs = 200\n",
    "x = a[b'Ekine']\n",
    "fig1, ax1 = plt.subplots()\n",
    "n, bins, patches = ax1.hist(x, nbs, density=True, alpha=0.75, label='E')\n",
    "ax1.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "cell_id": "00005-cf0c295b-5872-44c2-8a4d-a674bbb1fb0a",
    "execution_millis": 181,
    "execution_start": 1605009816015,
    "output_cleared": false,
    "source_hash": "4dcc46f7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of elements:  109619\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d4a5c4a459d4b9aab166ad4673c1ca6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc6e2366dc0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot the X,Y \n",
    "nbs = 200\n",
    "x = a[b'X']\n",
    "y = a[b'Y']\n",
    "print('Number of elements: ', len(x))\n",
    "# only keep 500 first elements \n",
    "n = 500\n",
    "x = x[:n]\n",
    "y = y[:n]\n",
    "fig2, ax2 = plt.subplots()\n",
    "ax2.scatter(x, y, alpha=0.75, label='x y')\n",
    "ax2.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "cell_id": "00006-caf99828-1f48-44af-92dc-0d08b4fdb2dc",
    "execution_millis": 1129,
    "execution_start": 1605009801896,
    "output_cleared": false,
    "source_hash": "51c65d4c"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c0ba76b793694a3fa24830aabb4bda90",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc6e22508b0>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot the theta angle (main direction dZ\n",
    "nbs = 500\n",
    "x = np.rad2deg(np.arccos(a[b'dZ']))\n",
    "fig3, ax3 = plt.subplots()\n",
    "n, bins, patches = ax3.hist(x, nbs, density=True, alpha=0.75, label='theta')\n",
    "ax3.set_xlim(160,180)\n",
    "ax3.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "cell_id": "00007-9738eca9-ea76-4f04-9949-5ebdba81fec6",
    "execution_millis": 548,
    "execution_start": 1605009803028,
    "output_cleared": false,
    "source_hash": "f91cf92a"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bc02c4bf169a40b49bbfbc9984439a17",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc6e3e61610>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot the phi angle\n",
    "nbs = 200\n",
    "x = np.rad2deg(np.arctan2(a[b'dY'], a[b'dX']))\n",
    "fig4, ax4 = plt.subplots()\n",
    "n, bins, patches = ax4.hist(x, nbs, density=True, alpha=0.75, label='phi')\n",
    "ax4.legend()"
   ]
  }
 ],
 "metadata": {
  "deepnote_execution_queue": [],
  "deepnote_notebook_id": "cffde93f-49bc-4a93-b8d3-146d7cad6d23",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
